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I. Motivation and Objective 

The recently developed second-order explicit and implicit total variation diminishing (TVD) shock- 
capturing methods of the Harten and Yee [1,2], Yee [3,4], and van Leer [5,6] types in conjunction 
with a generalized Roe’s approximate Riemann solver of Vinokur [7] and the generalized flux-vector 
splittings of Vinokur and Montagne [8] for two-dimensional hypersonic real gas flows are studied. A 
previous study [9] on one-dimensional unsteady problems indicated that these schemes produce good 
shock-capturing capability and that the state equation does not have a large effect on the general 
behavior of these methods for a wide range of flow conditions for equilibrium air. The objective of this 
paper is to investigate the applicability and shock resolution of these schemes for two-dimensional 
steady-state hypersonic blunt body flows. 

The main contribution of this paper is to identify some of the elements and parameters which can 
affect the convergence rate for high Mach numbers or real gases but have negligible effect for low 
Mach number cases for steady-state inviscid blunt body flows. In order to investigate these different 
points, two kinds of flows are considered. The blunt body calculations at Mach numbers higher than 
15 allow significant real gas effects to occur, while the case of an impinging shock provides a test on 
the treatment of slip surfaces and complex shock structures. In separate papers, a temporally second- 
order, implicit, time-accurate TVD-type algorithm for viscous steady and unsteady flows is studied. 
Studies show that the behavior of the schemes with various temporal differencing but similar spatial 
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discretization for inviscid and viscous flows are very different — in terms of stability and convergence 
rate. This point will be addressed in reference [10]. However, this paper only concerns itself with 
steady-state inviscid computations. 

In the following section, the generalized Roe’s approximate Riemann solver and flux- vector split- 
tings for real gases are reviewed. Due to space limitation, only the ADI linearized conservative 
implicit version of the Harten and Yee schemes [111 and Yee [3] is reviewed here since most of the il- 
lustrations are computed with this particular algorithm. The findings concerning the various aspects 
in improving the convergence rate and numerical examples are discussed in the subsequent sections. 

II. Description of the Numerical Algorithm 


The conservation laws for the two-dimensional Euler equations can be written in the form 


d£ dF(U) dG{U) 
dt dx ^ dy 


( 1 ) 


where U = [ p, m, n, e ] T , F = [ pu, mu + p, nu, eu+pu ] , and G = [ pv, mv, nv+p, ev-j-pv ] T . 
Here p is the density, m = pu is the x-component of the momentum per unit volume, n = pv is the 
y-component of the momentum per unit volume, p is the pressure, e = p[ e+ \ {u 2 + u 2 )] is the total 
internal energy per unit volume, and t is the specific internal energy. 


A generalized coordinate transformation of the form f = £(x,y) and p = p(x,y) which maintains 
the strong conservation-law form of equation (1) is given by 


dU_ dF(U) dd{U) 
dt ^ dp 


= 0 , 


( 2 ) 


where U = U / J , F = (£ X F + £ y G)/J, G = (p x F + p y G)/J, and J = £ x p y - £ y p x , the Jacobian 
transformation. Let A = dF/dU and B — dG/dU. Then the Jacobians A = dF/dU and B — 
dG / dU can be written as 

A = (£ X A + £ V B) (3a) 

B = (p x A + p y B). (3b) 


2.1. Riemann Solvers 

Here the usual approach of applying the one-dimensional scalar TVD schemes via the so called 
Riemann solvers for each direction in multidimensional nonlinear systems of hyperbolic conservation 
laws (see for example reference [2]) is used. The eigenvalues and eigenvectors of the Jacobian matrices 
A and B are used in approximate Riemann solvers. Given two states whose difference is At/, Roe 
[12] obtained an average A in the ^-direction, for example, satisfying A F = AAU for a perfect gas. 
The generalization by Vinokur [7] for an arbitrary gas involves the pressure derivatives x = ( dp/dp)~ 
and k = {dp/dt) , where t = pe. The relation c 2 = \ then gives the speed of sound, where 

h — e + p/ p. Introducing H = h + (u 2 + v 2 )/2, Vinokur found the same expressions for u, v and H 
as for the perfect gas, and that x and k must satisfy 

xAp + KAe = Ap. (4) 

Unique values of x and k are obtained by projecting the arithmetic averages of the values for the 
two states into this relation (see references [7] and [2] for the exact formulas). 
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Flux-vector splitting methods divide the flux F into several parts, each of which has a Jacobian 
matrix whose eigenvalues are all of one sign. The approach by Steger and Warming [13] made use 
of the relation F = AU , valid for a perfect gas. Van Leer [ 6 ] constructed a different splitting in 
which the eigenvalues of the split-flux Jacobians are continuous and one of them vanishes leading 
to sharper capture of transonic shocks. Vinokur and Montagne [ 8 ] showed that the expressions for 
both these splittings can be generalized to an arbitrary gas by using the variable 7 = pc 2 / p, and 
adding to the split energy flux a term equal to the product of the split mass flux and the quantity 
£ - C 2 j [ 7(7 - 1)] (see references [ 8 ] and [ 2 ] for the exact formulas). 


2.2 Description of the Implicit TVD schemes 

Let At be the time step and let the grid spacing be denoted by Af and Ar) such that £ = jAf 
and p = kAp. An implicit second-order in space, first-order in time TVD algorithm in generalized 
coordinates of Yee and Harten for two-dimensional systems (1) [2-4] can be written as 


frn+l , At r£ n+ i pr»+ 1 1 

u 3 a + Xe [ j+m 


+ 


At 

Ap 


\G n 

7 , 


n-fl 




(5) 


The functions F J+ 1 k and G hk+ i are the numerical fluxes in the £- and ^-directions evaluated at 
(j + h k) and [j, k + |), respectively. Typically, k can be expressed as 



~{F h k + F j+ i tk + 


( 6 ) 


Here i is the eigenvector matrix for dF/dU evaluated at some symmetric average of V and 
EW (for example, Roe average [12] for a perfect gas and generalized Roe average of Vinokur [7] 
for real gases). Similarly, one can define the numerical flux 1 in this manner. For viscous steady 
and unsteady flows, a fully implicit second-order in time and space algorithm (with the same spatial 
differencing for the convection terms) appears to be more stable and efficient (in terms of convergence 
rate) than (5). See references [10,14,15] for details. 


Second-order Symmetric TVD Scheme: The elements of the <J> J+ i in the ^-direction denoted by 

(<p l i ) for a spatially second-order symmetric TVD scheme [3,4] are 

J t 2 





(7a) 


The value a^ +i is the characteristic speed a 1 for dF/dU evaluated at some average between Uj ^ 
and Uj + i t k ■ The function ip is 


ip{z) 


\ z \ l 2 l ^ 

(z 2 + < 5 ! 2 )/ 25 i \z\<S 1 - 


(7b) 


Here t p(z) in equation (7b) is an entropy correction to |z| where is a small positive parameter. 
For steady-state problems containing strong shock waves, a proper control of the size of is very 
important, especially for hypersonic blunt-body flows. See reference [2] or section III for a discussion. 

An example of limiter function Q l 2 used in calculations is: 

2 



= minmod 


a 1 x ,a l 
3~i ’ 3- 



(7 c) 
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The minmod function of a list of arguments is equal to the smallest number in absolute value if the 
list of arguments is of the same sign, or is equal to zero if any arguments are of opposite sign. Here 

a 1 , are elements of 

J+2 


+ j Rj + l(Uj + l,lc Uj,k)- 


( 8 ) 


Second-Order Upwind TVD Scheme: The elements of the <$ J+ i in the ^-direction denoted by 


[(j> 1 i ) for a spatially second-order upwind TVD scheme [11,2] are 

J+ 2 


1 


K+i) = 2 ^K+ *)(»!■+» + ti) - + t£+ 4 K 


where 




iV'K+i) 


(»>+ 1 - ® i )/«*- +4 


+ 7 i+ 4 K+ *• 

(9a) 

<4 +i = 0 ' 

(9b) 


An example of limiter function g l 3 used in calculations is 


g\ — minmod a' _ , , a' , 

^ J 2 * ' 2 


(9c) 


^4 Conservative Linearized ADI Form for Steady-State Applications: A conservative linearized ADI 
form of equation (5) used mainly for steady-state applications as described in detail in references 
[3,11], can be written as 


J , iri _ 

/+ A£ * + Af i-4.* 


E" = - 


At 

A£l 


77»n 77’ n 


At 

Ar, 


/’* Ti /-vn 

J,fc+ 4 G J,k~i 


where 


The nonstandard notation 


/+ k+l - ~H\ , 

Arj JA+ 2 2 

C/ rt+1 = U n + E, 


E = E* , 


EJS 

2 


2 [^4j+x,fc + > 


H n 


l k+i = M, k+1 -w 






is used, and , , , , D' 7 , , , can be taken as 
3 + 2» fc 


(10a) 

(10b) 

(10c) 

(10d) 

(lOe) 

(lOf) 


n y+ 4 .*^* = i? J + 4 A di ag[V’(a' + r 

n i.fc+ 4 £7= *J.fc+i dia sM a U 4 


)]^ J ,fc + i(^V. fc + 1 



(10 B ) 

(lOh) 


Here A ;+ i,fc and fc+ i are (3) evaluated at [j + l,k) and (jf, k + l), respectively. The nonconserver- 
ative linearized implicit form suitable for steady-state calculations [2] is also considered. Numerical 
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study indicated that the latter form appears to be slightly less efficient in terms of convergence rate 
than the linearized conservative form. 


III. Enhancement of Convergence Rate for Hypersonic Flows 

The current study indicated that the following three elements can affect the convergence rate 
at hypersonic speeds: (a) the choice of the entropy correction parameter 6\, (b) the choice of the 
dependent variables on which the limiters are applied, and (c) the prevention of unphysical solutions 
during the initial transient stage. 

(a) . For blunt-body steady-state flows with M > 4, the initial flow conditions at the wall are 
obtained using the known wall temperature in conjunction with pressures computed from a modified 
Newtonian expression. Also, for implicit methods a slow startup procedure from freestream boundary 
conditions is necessary. Most importantly, it is advisable to use in equation (7b) as a function of 
the velocity and sound speed. In particular 

(£i) j + a = £(|u J+ i| + |v,- + i| + c y+4 ) (11a) 

(<$i) fc+ i = *(! u fc+il + K+il + c k+i) ( llb ) 

with 0.05 < 6 < 0.25 appear to be sufficient for the blunt-body flows for 4 < M < 25. Equation 
(11) is written in Cartesian coordinates. In the case of generalized coordinates, the u and t> should 
be replaced by the contravariant velocity components, and one half of the sound speed would be 
from the ^-direction and the other half would be from the ^-direction. For implicit methods, it is 
very important to use (11) in \p(z) on both the implicit and explicit operators. For the implicit 
operator, numerical experiments showed that the linearized conservative form (10) converges slightly 
faster than the linearized nonconservative form [11]. It seems also that when the freestream Mach 
number increases, the convergence rate of the linerarized conservative form (10) is slightly better 
than a simplified version which replaces Q^ +1 k and 0^ fc+i of (10g,h) by max; V’fay+i) and 
maxi V'( a / t+ x ) times the identity matrix. 

(b) . Higher-order TVD schemes in general involve limiter functions. However, there are options 
in choosing the types of dependent variables in applying limiters for system cases, in particular 
for systems in generalized coordinates. The choice of the dependent variables on which limiters are 
applied can affect the convergence process. In particular, due to the nonuniqueness of the eigenvectors 
R J+ 1 , the choice of the characteristic variables on which the limiters are applied play an important 
role in the convergence rate as the Mach number increases. For moderate Mach numbers, the different 
choice of the eigenvectors have negligible affect on the convergence rate. However, for large Mach 
number cases, the magnitudes of all the variables at the jump of the bow shock are not the same. In 
general, the jumps are much larger for the pressures than for the densities or total energy. Studies 
indicated that employing the form Rj+i such that the variation of the a are of the same order of 
magnitude as for the pressure would be a good choice for hypersonic flows. The form similar to the 
one used by Gnoffo [16] or Roe and Pike [17] can improve the convergence rate over the ones used 
in references [4,18]. 

(c) . Due to the large gradients and to the fact that the initial conditions are far from the steady- 
state physical solution, the path used by the implicit method can go through states with negative 
pressures if a large time step is employed. A convenient way to overcome the difficulties is to fix a 
minimum allowed value for the density and the pressure. With this safety check, the scheme allows 
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a much larger time step and converges several times faster. In addition, since the Roe’s average 
state allows the square of the average sound speed c 1 , to lie outside the interval between c 2 and 

3 ' 2 

c 2 +1 , C j + i might be negative even though c 2 and c 2 +1 are positive during the transient stage when 

the initial conditions are far from the steady-state physical solution. In this case, we replace c 2 , 

2 

by max(c 2 .,min(c 2 ,c 2 +1 )). This later safety check is in particular helpful for the symmetric TVD 
algorithm (7). 


IV. Numerical Results 

The current study on the shock resolution of the various schemes for two-dimensional steady-state 
blunt-body computations indicates similar trends as the one dimensional study [9]. The main issue 
appears to be their relative efficiency. Due to extra evaluations per dimension in the curve fitting 
between the left and right states in a real gas for the van Leer formulation, additional computation 
is required for the van Leer type schemes than the Harten and Yee [1,2], and Yee [3,4] types of TVD 
schemes. Here van Leer type schemes refer to the use of the MUSCL approach in conjunction with 
Roe type approximate Riemann solver [12] or flux-vector splittings [6,13]. Moreover, for steady-state 
applications, implicit methods are preferred over explicit methods because of the faster convergence 
rate. In addition, it is easier to obtain a noniterative linearized implicit operator for the Harten and 
Yee, and Yee type schemes than for the van Leer type schemes. For these reasons, the linearized 
implicit versions of Harten and Yee [ 1 1 ] and Yee [3] are preferred over the van Leer type schemes. 

Resolution of First- and Second-Order schemes : For problems containing complex shock structures, 
first-order upwind TVD schemes are too diffusive unless extremely fine grids are used. For a blunt- 
body flow containing a single steady bow shock only, the shock-capturing capability of a first-order 
upwind TVD scheme seems to be quite adequate if one is interested in the shock resolution only. 
However, a careful examination of the overall flow field of the density and Mach number contours of 
the first- and second-order TVD schemes compared with the exact solution (shock-fitting solution) 
reveals the inaccuracy of the first-order scheme. Figure 1 compares the resolution of the first-order 
(setting g l - = 0) and second-order upwind TVD schemes (10) using the Roe approximate Riemann 
solver [12] with the “exact solution” for a perfect gas (7 = 1.4) at a freestream Mach number of 10. 
The computations are performed on a 61 X 33 adapted grid for the full (half) cylinder, which yields 
a fairly good bow shock resolution by both schemes. However, the contour levels near the body are 
significantly shifted with the first-order scheme, while the second-order scheme reproduces almost 
identical results as the exact solution. 

Convergence Rate of Explicit and Implicit TVD Schemes at Hypersonic Speed: The five dif- 
ferent second-order TVD methods previously studied [9] in one dimension yield very similar shock- 
resolution for the blunt-body problem. In particular, for an inviscid blunt-body flow in the hypersonic 
equilibrium real gas range, the explicit second-order Harten and Yee, and Yee-Roe-Davis type TVD 
schemes [2-4] using the generalized approximate Riemann solver [7] produce similar shock-resolution 
but converge slightly faster than an explicit second-order van Leer type scheme using the generalized 
van Leer flux- vector splitting [8]. 

The freestream conditions for the current study are = 15 and 25, p ^ = 1.22 x 10 3 N/m 2 , 
Poo = 1.88 -2 kg/m 3 , and Too = 226°K. The grid size is 61 x 33 for the full (half) cylinder (figure 2). 
For the = 25 case, the shock stand off distance is at approximately fourteen points from the wall 
on the symmetry axis. The relaxation procedure for the explicit methods employs a second-order 
Runge-Kutta time-discretization with a CFL of 0.5 (solution not shown). The parameter 6 is set to 
a constant value of 0.15. Pressure and Mach number contours converge and stabilize after 3000-4000 
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steps but the convergence rate is much slower for the density (with a 2-3 order of magnitude drop 
in L2-norm residual). The bow shock is captured in two to three grid points. The curve fits of 
Srinivasan et al. [19] are used to generate the thermodynamic properties of the gas. 

The same flow condition was tested on the implicit scheme (10). The convergence rate is many 
times faster. Figures (3) and (4) show the Mach number, density, pressure and k contours computed 
by the linearized conservative ADI form of the upwind scheme (10) for Mach numbers 15 and 25. 
Figure 5 shows the slight advantage of the convergence rate of the linearized conservative implicit 
TVD scheme (10) over the linearized nonconservative implicit TVD scheme suggested in reference 
[11]. The convergence rate and shock resolution for the symmetric TVD scheme (10,7) behave 
similarly. For M <*, = 15 case, the L2-norm residual stagnated after a drop of four orders of magnitude. 
In general, for a perfect gas with 10 < < 25 and a not highly clustered grid, steady-state solutions 

can be reached in 800 steps with 12 orders of magnitude drop in the L 2 -norm residual. However, the 
convergence rate is at least twice as slow for the real gas counter-part. An important observation 
for the behavior of the convergence rate for the Mach 15 real gas case is that the discontinuities of 
the thermodynamic derivatives which exist in the curve fits of Srinivasan et al. [19] might be the 
major contributing factor. This is evident from figures (3d) and (4d) and from comparing with the 
convergence rate for the perfect gas result. 

Computations of impinging shocks: Figure (6) shows the Mach contours computed by the implicit 
upwind TVD scheme (10) of an inviscid shock-on-shock interaction on a blunt body in the low 
hypersonic range. Extensive study on flow fields of this type were reported in references [20-22] 
for the viscous case. This flow field is typical of what will be experienced by the inlet cowl of the 
National Aerospace Plane (NASP). The freestream conditions for this flow field are M a 0 = 4.6, 
= 14.93 N/m 2 , = 167°K, T w — 556°K, and 7 = 1.4 for a perfect gas. An oblique shock with 

an angle of 20.9° relative to the free stream impinges on the bow shock. Various types of interactions 
occur depending on where the impingement point is located on the bow shock. As shown by the 
Mach contours, the impinging shock has caused the stagnation point to be moved away from its 
undisturbed location at the symmetry line. The surface pressures at the new stagnation point can 
be several times larger than those at the undisturbed location of the stagnation point. In addition, a 
slip surface emanates from the bow shock and impinging shock intersection point and is intercepted 
by a shock wave which starts at the upper kink of the bow shock. The interacting shock waves and 
slip surfaces are confined to a very small region and must be captured accurately by the numerical 
scheme if the proper surface pressures and heat transfer rates are to be predicted correctly. The 
77 x 77 grid used and the convergence rate computed by the implicit scheme (10) are shown in figure 
(6). Though the pattern of the flow is significantly more complicated than for the previous cases, the 
convergence rate remains quite satisfactory. Detailed study of viscous steady and unsteady flow fields 
of this type using a fully implicit second-order time-accurate scheme [10] of the same numerical flux 
(6-9) for the convection terms are reported in [10,14,15]. It was found that for viscous computations, 
the scheme suggested by Yee et al. [10] is more robust than equation (10) which is best suited for 
steady-state inviscid flows. 


IV. Concluding Remarks 

Some numerical aspects of the TVD schemes that can affect the convergence rate for hypersonic 
Mach numbers or real gas flows but have negligible effect on low Mach number or perfect gas flows 
are identified. Improvements have been made to the various TVD algorithms to speed up the 
convergence rate in the hypersonic flow regime. Even with the improvement though, the convergence 
is in general slightly slower for a real gas than for a perfect gas. The nonsmoothness in the curve fits 
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of Srinivasan et al. may be a major contributing factor in slowing down the convergence rate. Due 
to extra evaluations per dimension in the curve fitting between the left and right states in a real gas 
for the van Leer formulation, more computation is required for the van Leer type schemes than for 
the Harten and Yee, and Yee types of TVD schemes. 

Aside from the difference in convergence rate, the numerical results confirm the findings of the 
one dimensional study. The different methods yield very similar shock-resolution on the blunt body 
problem with freestream Mach numbers up to 25, and the state equation does not have a large effect 
on the general behavior of these methods. Further improvements on the ADI relaxation algorithm 
could speed up the convergence rate even more. 
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Fig. 3 The Mach contours (a), density contours (b), pressure contours (c) and k (d) computed by 
second-order implicit TVD scheme for an equilibrium real gas at M a 0 = 15. 
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The Mach contours (a), density contours (b), pressure contours (c) and k (d) computed by a 
second-order implicit TVD scheme for an equilibrium real gas at = 25. 
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Fig. 5 Comparison of the Z, 2 -norm residual of a linearized conservative implicit operator (a) and a 
linearized nonconservative implicit operator (b) for an equilibrium real gas at M <*, = 25. 
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